//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2009 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : various
// Module     : ????
// File       : qshep2d.cpp
//
// NOSTYLE    : this file is an only very slightly upgraded output of f2c !
// TODO       : make this from plain and ugly C into nice C++ code sometime? ASR
// Original   : qshep2d.f90 (http://people.sc.fsu.edu/~burkardt/f_src/qshep2d/qshep2d.html)
//----------------------------------------------------------------------------

#include <math.h>
#include "qshep2d_c.h"

#define min(a, b) ((a)<(b)?(a):(b))
#define max(a, b) ((a)>(b)?(a):(b))
#define abs(a) ((a)<0?-(a):(a))

/* *****************************************************************************80 */

/* ! GETNP2 seeks the closest unmarked node to a point. */

/*  Discussion: */

/*    GETNP2 uses the cell method to find the closest unmarked node NP */
/*    to a specified point P, given a set of N nodes and the data structure */
/*    defined by STORE2. */

/*    NP is then marked by negating LNEXT(NP).  Thus, the closest M nodes to */
/*    P may be determined by a sequence of M calls to this routine. */

/*    If the point P is itself actually a node K, and you want to find the */
/*    nearest point to P that is not node K, then you must be sure to mark */
/*    node K before calling. */

/*    The search is begun in the cell containing or closest to P and proceeds */
/*    outward in rectangular layers until all cells which contain points */
/*    within distance R of P have been searched.  R is the distance from P to */
/*    the first unmarked node encountered, or infinite if no unmarked nodes */
/*    are present. */

/*    Input parameters other than LNEXT are not altered by this routine. */
/*    With the exception of ( PX, PY ) and the signs of LNEXT elements, */
/*    these parameters should be unaltered from their values on output */
/*    from subroutine STORE2. */

/*  Modified: */

/*    10 July 1999 */

/*  Author: */

/*    Robert Renka, */
/*    University of North Texas */

/*  Reference: */

/*    Robert Renka, */
/*    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate */
/*    interpolation of scattered data, */
/*    ACM Transactions on Mathematical Software, */
/*    Volume 14, 1988, pages 149-150. */

/*  Parameters: */

/*    Input, double precision PX, PY, the (X,Y) coordinates of the point P */
/*    whose nearest unmarked neighbor is to be found. */

/*    Input, double precision X(N), Y(N), the coordinates of the nodes at which */
/*    data has been supplied. */

/*    Input, int NR, the number of rows and columns in the cell grid. */
/*    NR must be at least 1. */

/*    Input, int LCELL(NR,NR), array of nodal indices associated */
/*    with cells. */

/*    Input/output, int LNEXT(N), contains next-node indices ( or their */
/*    negatives ).  On return, if the output value of NP is nonzero, then */
/*    LNEXT(NP) will be negative. */

/*    Input, double precision XMIN, YMIN, DX, DY, the minimum nodal X, Y */
/*    coordinates, and the X, Y dimensions of a cell.  DX and DY must be */
/*    positive. */

/*    Output, int NP, the index into the vectors X and Y of the nearest */
/*    unmarked node to the point P.  NP will be 0 if all nodes are marked */
/*    or if the values of NR, DX, DY are illegal.  LNEXT(NP) will be less */
/*    than 0 if NP is nonzero (this marks node NP as being used now). */

/*    Output, double precision DSQ, if NP is nonzero, then DSQ is the */
/*    squared distance between P and node NP. */

/*  Local Parameters: */

/*    first = true iff the first unmarked node has yet to be encountered, */

/*    imin, imax, jmin, jmax = cell indices defining the range of the search, */

/*    delx, dely = px-xmin and py-ymin, */

/*    i0, j0 = cell containing or closest to P, */

/*    i1, i2, j1, j2 = cell indices of the layer whose intersection with */

/*    the range defined by imin,...,jmax is currently being searched. */

static int getnp2(double *px, double *py, double *x,
                  double *y, int *nr, int *lcell, int *lnext,
                  double *xmin, double *ymin, double *dx, double *dy,
                  int *np, double *dsq)
{
    /* System generated locals */
    int lcell_dim1, lcell_offset, i__1, i__2;
    double d__1, d__2;

    /* Local variables */
    static int i__, j, l;
    static double r__;
    static int i0, i1, i2, j0, j1, j2, ln;
    static double xp, yp, rsq, delx, dely;
    static int imax, imin, jmax, jmin, lmin;
    static bool first;
    static double rsmin;

    /* Parameter adjustments */
    --x;
    --y;
    lcell_dim1 = *nr;
    lcell_offset = 1 + lcell_dim1;
    lcell -= lcell_offset;
    --lnext;

    /* Function Body */
    xp = *px;
    yp = *py;

    /*  Test for invalid input parameters. */

    if (*nr < 1 || *dx <= 0. || *dy <= 0.)
    {
        *np = 0;
        *dsq = 0.;
    }

    /*  Initialize parameters: */

    first = true;
    imin = 1;
    imax = *nr;
    jmin = 1;
    jmax = *nr;
    delx = xp - *xmin;
    dely = yp - *ymin;
    i0 = (int)(delx / *dx) + 1;
    i0 = max(i0, 1);
    i0 = min(i0, *nr);
    j0 = (int)(dely / *dy) + 1;
    j0 = max(j0, 1);
    j0 = min(j0, *nr);
    i1 = i0;
    i2 = i0;
    j1 = j0;
    j2 = j0;

    /*  Outer loop on layers, inner loop on layer cells, excluding */
    /*  those outside the range (imin,imax) x (jmin,jmax). */

L1:
    i__1 = j2;
    for (j = j1; j <= i__1; ++j)
    {
        if (jmax < j)
        {
            /*      exit */
            goto L201;
        }
        if (j < jmin)
        {
            /* cycle */
            goto L301;
        }
        i__2 = i2;
        for (i__ = i1; i__ <= i__2; ++i__)
        {
            if (imax < i__)
            {
                /* exit */
                goto L202;
            }
            if (i__ < imin)
            {
                /*        cycle */
                goto L302;
            }
            if (j != j1 && j != j2 && i__ != i1 && i__ != i2)
            {
                /*        cycle */
                goto L302;
            }

            /*  Search cell (i,j) for unmarked nodes l. */

            l = lcell[i__ + j * lcell_dim1];
            if (0 < l)
            {

                /*  Loop on nodes in cell (i,j). */

L2:
                ln = lnext[l];

                /*  Node L is the first unmarked neighbor of P encountered. */

                /*  Initialize LMIN to the current candidate for np, and */
                /*  rsmin to the squared distance from p to lmin.  imin, */
                /*  imax, jmin, and jmax are updated to define the smal- */
                /*  lest rectangle containing a circle of radius r = */
                /*  sqrt(rsmin) centered at p, and contained in (1,nr) x */
                /*  (1,nr) (except that, if p is outside the rectangle */
                /*  defined by the nodes, it is possible that imin .gt. */
                /*  nr, imax < 1, nr < jmin, or jmax < 1). */

                if (0 <= ln)
                {
                    /* Computing 2nd power */
                    d__1 = x[l] - xp;
                    /* Computing 2nd power */
                    d__2 = y[l] - yp;
                    rsq = d__1 * d__1 + d__2 * d__2;
                    if (first)
                    {
                        lmin = l;
                        rsmin = rsq;
                        r__ = sqrt(rsmin);
                        imin = (int)((delx - r__) / *dx) + 1;
                        imin = max(imin, 1);
                        imax = (int)((delx + r__) / *dx) + 1;
                        imax = min(imax, *nr);
                        jmin = (int)((dely - r__) / *dy) + 1;
                        jmin = max(jmin, 1);
                        jmax = (int)((dely + r__) / *dy) + 1;
                        jmax = min(jmax, *nr);
                        first = false;
                    }
                    else
                    {
                        if (rsq < rsmin)
                        {
                            lmin = l;
                            rsmin = rsq;
                        }
                    }
                }
                if (abs(ln) != l)
                {
                    l = abs(ln);
                    goto L2;
                }
            }
            /*   end do */
L302:
            /* L102: */
            ;
        }
L202:
        /*  end do */
        L301 :
            /* L101: */
            ;
    }
L201:

    /*  Test for termination of loop on cell layers. */

    if (imin < i1 || i2 < imax || jmin < j1 || j2 < jmax)
    {
        --i1;
        ++i2;
        --j1;
        ++j2;
        goto L1;
    }
    if (first)
    {
        *np = 0;
        *dsq = 0.;
    }
    else
    {
        *np = lmin;
        *dsq = rsmin;
        lnext[lmin] = -lnext[lmin];
    }
    return 0;
} /* getnp2_ */


/* *****************************************************************************80 */

/* ! ROTATE applies a Givens rotation. */

/*  Discussion: */

/*    The rotation has the form: */

/*      (   C  S ) */
/*      ( - S  C ) */

/*    and is essentially applied to a 2 by N matrix: */

/*      ( X(1) X(2) ... X(N) ) */
/*      ( Y(1) Y(2) ... Y(N) ) */

/*  Modified: */

/*    28 June 1999 */

/*  Parameters: */

/*    Input, int N, the dimension of the vectors. */

/*    Input, double precision C, S, the cosine and sine entries of the Givens */
/*    rotation matrix.  These may be determined by subroutine GIVENS. */

/*    Input/output, double precision X(N), Y(N), the rotated vectors. */

static void rotate(int *n, double *c__, double *s, double *x, double *y)
{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int i__;
    static double xi, yi;

    /* Parameter adjustments */
    --y;
    --x;

    /* Function Body */
    if (*n <= 0)
    {
        return;
    }
    else if (*c__ == 1. && *s == 0.)
    {
        return;
    }
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        xi = x[i__];
        yi = y[i__];
        x[i__] = *c__ * xi + *s * yi;
        y[i__] = -(*s) * xi + *c__ * yi;
        /*  end do */
        /* L115: */
    }
} /* rotate_ */


/* *****************************************************************************80 */

/* ! SETUP2 sets up a row of the least squares regression matrix. */

/*  Discussion: */

/*    SETUP2 sets up the I-th row of an augmented regression matrix for */
/*    a weighted least-squares fit of a quadratic function Q(X,Y) to a set */
/*    of data values F, where Q(XK,YK) = FK. */

/*    The first 3 columns are quadratic terms, and are scaled by 1/S2. */
/*    The fourth and fifth columns represent linear terms, and are scaled */
/*    by 1/S1. */

/*    If D = 0, or R <= D, the weight is */
/*      0, */
/*    else if D < R, the weight is */
/*      (R-D)/(R*D), */
/*    where D is the distance between nodes I and K, and R is a maximum */
/*    influence distance. */

/*  Modified: */

/*    05 July 1999 */

/*  Author: */

/*    Robert Renka, */
/*    University of North Texas */

/*  Reference: */

/*    Robert Renka, */
/*    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate */
/*    interpolation of scattered data, */
/*    ACM Transactions on Mathematical Software, */
/*    Volume 14, 1988, pages 149-150. */

/*  Parameters: */

/*    Input, double precision XK, YK, FK, the coordinates and value of the data */
/*    at data node K. */

/*    Input, double precision XI, YI, FI, the coorindates and value of the data */
/*    at data node I. */

/*    Input, double precision S1, S2, reciprocals of the scale factors. */

/*    Input, double precision R, the maximum radius of influence about node K. */

/*    Output, double precision ROW(6), a row of the augmented regression matrix. */

static void setup2(double *xk, double *yk, double *fk,
                   double *xi, double *yi, double *fi, double *s1,
                   double *s2, double *r__, double *row)
{
    /* Local variables */
    static double d__, w, dx, dy;
    static int ijk;

    /* Parameter adjustments */
    --row;

    /* Function Body */
    dx = *xi - *xk;
    dy = *yi - *yk;
    d__ = sqrt(dx * dx + dy * dy);
    if (d__ <= 0. || *r__ <= d__)
    {
        for (ijk = 1; ijk <= 6; ++ijk)
        {
            row[ijk] = 0.;
            /* L100: */
        }
        /*    row(1:6) = 0.0D+00 */
    }
    else
    {
        w = (*r__ - d__) / *r__ / d__;
        row[1] = dx * dx * w / *s2;
        row[2] = dx * dy * w / *s2;
        row[3] = dy * dy * w / *s2;
        row[4] = dx * w / *s1;
        row[5] = dy * w / *s1;
        row[6] = (*fi - *fk) * w;
    }
} /* setup2_ */


/* *****************************************************************************80 */

/* ! STORE2 creates a cell data structure for the scattered data. */

/*  Discussion: */

/*    STORE2 is given a set of N arbitrarily distributed nodes in the */
/*    plane and creates a data structure for a cell-based method of */
/*    solving closest-point problems.  The smallest rectangle containing */
/*    all the nodes is partitioned into an NR by NR uniform grid of cells, */
/*    and nodes are associated with cells. */

/*    In particular, the data structure stores the indices of the nodes */
/*    contained in each cell.  For a uniform random distribution of nodes, */
/*    the nearest node to an arbitrary point can be determined in constant */
/*    expected time. */

/*  Modified: */

/*    05 July 1999 */

/*  Author: */

/*    Robert Renka */
/*    University of North Texas */

/*  Reference: */

/*    Robert Renka, */
/*    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate */
/*    interpolation of scattered data, */
/*    ACM Transactions on Mathematical Software, */
/*    Volume 14, 1988, pages 149-150. */

/*  Parameters: */

/*    Input, int N, the number of data nodes.  N must be at least 2. */

/*    Input, double precision X(N), Y(N), the coordinates of the data nodes. */

/*    Input, int NR, the number of rows and columns in the grid.  The */
/*    cell density, or average number of data nodes per cell, is */
/*      D = N / ( NR * NR ). */
/*    A recommended value, based on empirical evidence, is */
/*      D = 3. */
/*    Hence, the corresponding value of NR is recommended to be about */
/*      NR = SQRT ( N / 3 ). */
/*    NR must be at least 1. */

/*    Output, int LCELL(NR,NR), an array set up so that LCELL(I,J) */
/*    contains the index (for X and Y) of the first data node (that is, the */
/*    data node with smallest index) in the (I,J) cell.  LCELL(I,J) will be 0 if */
/*    no data nodes are contained in the (I,J) cell.  The upper right corner of */
/*    the (I,J) cell has coordinates */
/*      ( XMIN + I * DX, YMIN + J * DY ). */

/*    Output, int LNEXT(N), an array of next-node indices.  LNEXT(K) */
/*    contains the index of the next node in the cell which contains node K, */
/*    or LNEXT(K) = K if K is the last node in the cell. */
/*    The data nodes contained in a cell are ordered by their indices. */
/*    If, for example, cell (I,J) contains nodes 2, 3, and 5 and no others, */
/*    then: */

/*      LCELL(I,J) = 2, (index of the first data node) */

/*      LNEXT(2) = 3, */
/*      LNEXT(3) = 5, */
/*      LNEXT(5) = 5. */

/*    Output, double precision XMIN, YMIN, the X, Y coordinates of the lower */
/*    left corner of the rectangle defined by the data nodes.  The upper right */
/*    corner is ( XMAX, YMAX ), where */
/*      XMAX = XMIN + NR * DX, */
/*      YMAX = YMIN + NR * DY. */

/*    Output, double precision DX, DY, the X and Y dimensions of the */
/*    individual cells. */
/*      DX = ( XMAX - XMIN ) / NR */
/*      DY = ( YMAX - YMIN ) / NR, */
/*    where XMIN, XMAX, YMIN and YMAX are the extrema of X and Y. */

/*    Output, int IER, an error indicator. */
/*    0, if no errors were encountered. */
/*    1, if N < 2 or NR < 1. */

/*    2, if DX = 0 or DY = 0. */

static void store2(int *n, double *x, double *y,
                   int *nr, int *lcell, int *lnext, double *xmin,
                   double *ymin, double *dx, double *dy, int *ier)
{
    /* System generated locals */
    int lcell_dim1, lcell_offset, i__1, i__2;

    /* Local variables */
    static int i__, j, k, l, ijk, kji;
    static double xmax, ymax;

    /* Parameter adjustments */
    --lnext;
    --y;
    --x;
    lcell_dim1 = *nr;
    lcell_offset = 1 + lcell_dim1;
    lcell -= lcell_offset;

    /* Function Body */
    *ier = 0;
    if (*n < 2)
    {
        *ier = 1;
        return;
    }
    if (*nr < 1)
    {
        *ier = 1;
        return;
    }

    /*  Compute the dimensions of the (X,Y) rectangle containing all the data nodes. */

    *xmin = x[1];
    xmax = x[1];
    *ymin = y[1];
    ymax = y[1];
    i__1 = *n;
    for (ijk = 1; ijk <= i__1; ++ijk)
    {
        if (x[ijk] < *xmin)
        {
            *xmin = x[ijk];
        }
        if (y[ijk] < *ymin)
        {
            *ymin = y[ijk];
        }
        if (x[ijk] > xmax)
        {
            xmax = x[ijk];
        }
        if (y[ijk] > ymax)
        {
            ymax = y[ijk];
        }
        /* L500: */
    }
    /*  xmin = minval ( x(1:n) ) */
    /*  xmax = maxval ( x(1:n) ) */
    /*  ymin = minval ( y(1:n) ) */
    /*  ymax = maxval ( y(1:n) ) */

    /*  Compute the dimensions of a single cell. */

    *dx = (xmax - *xmin) / *nr;
    /* real ( nr, kind = 8 ) */
    *dy = (ymax - *ymin) / *nr;

    /*  Test for zero area. */

    /* real ( nr, kind = 8 ) */
    if (*dx == 0. || *dy == 0.)
    {
        *ier = 2;
        return;
    }

    /*  Initialize LCELL. */

    i__1 = *nr;
    for (ijk = 1; ijk <= i__1; ++ijk)
    {
        i__2 = *nr;
        for (kji = 1; kji <= i__2; ++kji)
        {
            lcell[ijk + kji * lcell_dim1] = 0;
            /* L601: */
        }
        /* L600: */
    }
    /*        lcell(1:nr,1:nr) = 0 */

    /*  Loop on nodes, storing indices in LCELL and LNEXT. */

    for (k = *n; k >= 1; --k)
    {
        i__ = (int)((x[k] - *xmin) / *dx) + 1;
        i__ = min(i__, *nr);
        j = (int)((y[k] - *ymin) / *dy) + 1;
        j = min(j, *nr);
        l = lcell[i__ + j * lcell_dim1];
        if (l != 0)
        {
            lnext[k] = l;
        }
        else
        {
            lnext[k] = k;
        }
        lcell[i__ + j * lcell_dim1] = k;
        /*  end do */
        /* L116: */
    }
} /* store2_ */


/* *****************************************************************************80 */

/* ! GIVENS constructs a Givens plane rotation. */

/*  Discussion: */

/*    The transformation has the form of a 2 by 2 matrix G(C,S): */

/*      (   C  S ) */
/*      ( - S  C ) */

/*    where C*C + S*S = 1, which zeroes the second entry of the */
/*    the column vector ( A, B ) when C and S are properly chosen. */
/*    A call to GIVENS is normally followed by a call to ROTATE */
/*    which computes the product of G(C,S) with a 2 by N matrix. */

/*  Modified: */

/*    10 July 1999 */

/*  Parameters: */

/*    Input/output, double precision A, B. */

/*    On input, A and B define the 2-vector whose second entry (B) is */
/*    to be annihilated by a Givens rotation. */

/*    On output, A has been overwritten by a value */
/*      R = +/- SQRT ( A*A + B*B ) */
/*    and B has been overwritten by a value Z which allows C */
/*    and S to be recovered as: */

/*      if | Z | <= 1, then */
/*        C = SQRT ( 1 - Z*Z ), */
/*        S = Z */
/*      else if 1 < | Z | then */
/*        C = 1 / Z, */
/*        S = SQRT ( 1 - C*C ). */

/*    Output, double precision C, S, the components of the Givens */
/*    transformation, which may be computed by: */

/*      C = +/- A / SQRT ( A*A + B*B ) */
/*      S = +/- B / SQRT ( A*A + B*B ) */

/*  Local parameters: */

/*  r =        c*a + s*b = +/-sqrt(a*a+b*b) */

/*  u,v =   variables used to scale a and b for computing r */

static void givens(double *a, double *b, double *c__, double *s)
{
    /* Local variables */
    static double r__, u, v;

    /*  abs(b) < abs(a) */

    /*  Note that r has the sign of a, 0 < c, and s has */
    /*  sign(a)*sign(b). */

    if (fabs(*b) < fabs(*a))
    {
        u = *a * 2.;
        v = *b / u;
        r__ = sqrt(v * v + .25) * u;
        *c__ = *a / r__;
        *s = v * 2. * *c__;
        *b = *s;
        *a = r__;

        /*  abs(a) <= abs(b) */

        /*  Store r in a. */
        /*  Note that r has the sign of b, 0 < s, and c has sign(a)*sign(b). */

    }
    else if (*b != 0.)
    {
        u = *b * 2.;
        v = *a / u;
        *a = sqrt(v * v + .25) * u;
        *s = *b / *a;
        *c__ = v * 2. * *s;
        if (*c__ != 0.)
        {
            *b = 1. / *c__;
        }
        else
        {
            *b = 1.;
        }

        /*  a = b = 0. */

    }
    else
    {
        *c__ = 1.;
        *s = 0.;
    }
} /* givens_ */


/* *****************************************************************************80 */

/* ! QSHEP2 computes an interpolant to scattered data in the plane. */

/*  Discussion: */

/*    QSHEP2 computes a set of parameters A and RSQ defining a smooth, */
/*    once continuously differentiable, bi-variate function Q(X,Y) which */
/*    interpolates given data values F at scattered nodes (X,Y). */

/*    The interpolant function Q(X,Y) may be evaluated at an arbitrary point */
/*    by passing the parameters A and RSQ to the function QS2VAL.  The */
/*    first derivatives dQdX(X,Y) and dQdY(X,Y) may be evaluated by */
/*    subroutine QS2GRD. */

/*    The interpolation scheme is a modified quadratic Shepard method: */

/*      Q = ( W(1) * Q(1) + W(2) * Q(2) + .. + W(N) * Q(N) ) */
/*        / ( W(1)        + W(2)        + .. + W(N) ) */

/*    for bivariate functions W(K) and Q(K).  The nodal functions are given by */

/*      Q(K)(X,Y) = */
/*          F(K) */
/*        + A(1,K) * ( X - X(K) )**2 */
/*        + A(2,K) * ( X - X(K) ) * ( Y - Y(K) ) */
/*        + A(3,K) * ( Y - Y(K) )**2 */
/*        + A(4,K) * ( X - X(K) ) */
/*        + A(5,K) * ( Y - Y(K) ). */

/*    Thus, Q(K) is a quadratic function which interpolates the */
/*    data value at node K.  Its coefficients A(*,K) are obtained */
/*    by a weighted least squares fit to the closest NQ data */
/*    points with weights similar to W(K).  Note that the radius */
/*    of influence for the least squares fit is fixed for each */
/*    K, but varies with K. */

/*    The weights are taken to be */

/*      W(K)(X,Y) = ( (R(K)-D(K))+ / R(K) * D(K) )**2 */

/*    where (R(K)-D(K))+ = 0 if R(K) <= D(K) and D(K)(X,Y) is */
/*    the euclidean distance between (X,Y) and (X(K),Y(K)).  The */
/*    radius of influence R(K) varies with K and is chosen so */
/*    that NW nodes are within the radius.  Note that W(K) is */
/*    not defined at node (X(K),Y(K)), but Q(X,Y) has limit F(K) */
/*    as (X,Y) approaches (X(K),Y(K)). */

/*  Author: */

/*    Robert Renka, */
/*    University of North Texas */

/*  Reference: */

/*    Robert Renka, */
/*    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate */
/*    interpolation of scattered data, */
/*    ACM Transactions on Mathematical Software, */
/*    Volume 14, 1988, pages 149-150. */

/*  Parameters: */

/*    Input, int N, the number of nodes (X,Y) at which data values */
/*    are given.  N must be at least 6. */

/*    Input, double precision X(N), Y(N), the coordinates of the nodes at which */
/*    data has been supplied. */

/*    Input, double precision F(N), the data values. */

/*    Input, int NQ, the number of data points to be used in the least */
/*    squares fit for coefficients defining the nodal functions Q(K). */
/*    A highly recommended value is NQ = 13. */
/*    NQ must be at least 5, and no greater than the minimum of 40 and N-1. */

/*    Input, int NW, the number of nodes within (and defining) the radii */
/*    of influence R(K) which enter into the weights W(K).  For N */
/*    sufficiently large, a recommended value is NW = 19.   NW must be */
/*    at least 1, and no greater than the minimum of 40 and N-1. */

/*    Input, int NR, the number of rows and columns in the cell grid */
/*    defined in subroutine STORE2.  A rectangle containing the nodes */
/*    is partitioned into cells in order to increase search efficiency. */
/*    NR = SQRT(N/3) is recommended.  NR must be at least 1. */

/*    Output, int LCELL(NR,NR), array of nodal indices associated */
/*    with cells. */

/*    Output, int LNEXT(N), contains next-node indices ( or their */
/*    negatives ). */

/*    Output, double precision XMIN, YMIN, DX, DY, the minimum nodal X, Y */
/*    coordinates, and the X, Y dimensions of a cell. */

/*    Output, double precision RMAX, the square root of the largest element */
/*    in RSQ, the maximum radius of influence. */

/*    Output, double precision RSQ(N), the squared radii which enter into */
/*    the weights defining the interpolant Q. */

/*    Output, double precision A(5,N), the coefficients for the nodal functions */
/*    defining the interpolant Q. */

/*    Output, int IER, error indicator. */
/*    0, if no errors were encountered. */
/*    1, if N, NQ, NW, or NR is out of range. */
/*    2, if duplicate nodes were encountered. */
/*    3, if all nodes are collinear. */

/*  Local parameters: */

/* av =        root-mean-square distance between k and the */
/*             nodes in the least squares system (unless */
/*             additional nodes are introduced for stabil- */
/*             ity).      the first 3 columns of the matrix */
/*             are scaled by 1/avsq, the last 2 by 1/av */
/* avsq =      av*av */
/* b =         transpose of the augmented regression matrix */
/* c =         first component of the plane rotation used to */
/*             zero the lower triangle of b**t -- computed */
/*             by subroutine givens */
/* ddx,ddy =   local variables for dx and dy */
/* dmin =      minimum of the magnitudes of the diagonal */
/*             elements of the regression matrix after */
/*             zeros are introduced below the diagonal */
/* DTOL =      tolerance for detecting an ill-conditioned */
/*             system.  the system is accepted when DTOL <= DMIN. */
/* fk =        data value at node k -- f(k) */
/* i =         index for a, b, and npts */
/* ib =        do-loop index for back solve */
/* ierr =      error flag for the call to store2 */
/* irow =      row index for b */
/* j =         index for a and b */
/* jp1 =       j+1 */
/* k =         nodal function index and column index for a */
/* lmax =      maximum number of npts elements (must be con- */
/*             sistent with the dimension statement above) */
/* lnp =       current length of npts */
/* neq =       number of equations in the least squares fit */
/* nn,nnr =    local copies of n and nr */
/* np =        npts element */
/* npts =      array containing the indices of a sequence of */
/*             nodes to be used in the least squares fit */
/*             or to compute rsq.  the nodes are ordered */
/*             by distance from k and the last element */
/*             (usually indexed by lnp) is used only to */
/*             determine rq, or rsq(k) if NQ < NW. */
/* nqwmax =    max(nq,nw) */
/* rq =        radius of influence which enters into the */
/*             weights for q(k) (see subroutine setup2) */
/* rs =        squared distance between k and npts(lnp) -- */
/*             used to compute rq and rsq(k) */
/* rsmx =      maximum rsq element encountered */
/* rsold =     squared distance between k and npts(lnp-1) -- */
/*             used to compute a relative change in rs */
/*             between succeeding npts elements */
/* RTOL =      tolerance for detecting a sufficiently large */
/*             relative change in rs.  if the change is */
/*             not greater than RTOL, the nodes are */
/*             treated as being the same distance from k */
/* rws =       current value of rsq(k) */
/* s =         second component of the plane givens rotation */
/* SF =        marquardt stabilization factor used to damp */
/*             out the first 3 solution components (second */
/*             partials of the quadratic) when the system */
/*             is ill-conditioned.  as SF increases, the */
/*             fitting function approaches a linear */
/* sum2 =      sum of squared euclidean distances between */
/*             node k and the nodes used in the least */
/*             squares fit (unless additional nodes are */
/*             added for stability) */
/* t =         temporary variable for accumulating a scalar */
/*             product in the back solve */
/* xk,yk =     coordinates of node k -- x(k), y(k) */

/* xmn,ymn =   local variables for xmin and ymin */

void qshep2(int n, double *x, double *y,
            double *f, int nq, int nw, int nr, int *lcell,
            int *lnext, double *xmin, double *ymin, double *dx,
            double *dy, double *rmax, double *rsq, double *a,
            int *ier)
{
    /* System generated locals */
    int lcell_dim1, lcell_offset, i__1, i__2, i__3;
    double d__1, d__2;

    /* Local variables */
    static double b[36] /* was [6][6] */, c__;
    static int i__, j, k;
    static double s, t, fk, av, sf;
    static int nn, np;
    static double rq, xk, rs, yk;
    static int jp1, ijk;
    static double ddx, ddy;
    static int neq, lnp, nnr;
    static double xmn, ymn, rws, sum2, dmin__;
    static int ierr;
    static double dtol;
    static int lmax;
    static double avsq;
    static int irow;
    static double rtol;
    static int npts[40];
    static double rsmx, rsold;
    static int nqwmax;


    /* Parameter adjustments */
    a -= 6;
    --rsq;
    --lnext;
    --f;
    --y;
    --x;
    lcell_dim1 = nr;
    lcell_offset = 1 + lcell_dim1;
    lcell -= lcell_offset;

    /* Function Body */
    if (n < 6)
    {
        *ier = 1;
        return;
    }
    dtol = .01;
    /* parameter */
    rtol = 1e-5;
    /* parameter */
    sf = 1.;
    /* parameter */
    nn = n;
    nnr = nr;
    nqwmax = max(nq, nw);
    /* Computing MIN */
    i__1 = 40, i__2 = n - 1;
    lmax = min(i__1, i__2);
    if (nq < 5)
    {
        *ier = 1;
        return;
    }
    else if (nw < 1)
    {
        *ier = 1;
        return;
    }
    else if (lmax < nqwmax)
    {
        *ier = 1;
        return;
    }
    else if (nr < 1)
    {
        *ier = 1;
        return;
    }

    /*  Create the cell data structure, and initialize RSMX. */

    store2(&nn, &x[1], &y[1], &nnr, &lcell[lcell_offset], &lnext[1], &xmn, &
           ymn, &ddx, &ddy, &ierr);
    if (ierr != 0)
    {
        *xmin = xmn;
        *ymin = ymn;
        *dx = ddx;
        *dy = ddy;
        *ier = 3;
        return;
    }
    rsmx = 0.;

    /*  Outer loop on node K. */

    i__1 = nn;
    for (k = 1; k <= i__1; ++k)
    {
        xk = x[k];
        yk = y[k];
        fk = f[k];

        /*  Mark node K to exclude it from the search for nearest neighbors. */

        lnext[k] = -lnext[k];

        /*  Initialize for loop on NPTS. */

        rs = 0.;
        sum2 = 0.;
        rws = 0.;
        rq = 0.;
        lnp = 0;

        /*  Compute NPTS, LNP, RWS, NEQ, RQ, and AVSQ. */

L1:
        sum2 += rs;
        if (lnp == lmax)
        {
            goto L3;
        }
        ++lnp;
        rsold = rs;
        getnp2(&xk, &yk, &x[1], &y[1], &nnr, &lcell[lcell_offset], &lnext[1],
               &xmn, &ymn, &ddx, &ddy, &np, &rs);
        if (rs == 0.)
        {
            *ier = 2;
            return;
        }
        npts[lnp - 1] = np;
        if ((rs - rsold) / rs < rtol)
        {
            goto L1;
        }
        if (rws == 0. && nw < lnp)
        {
            rws = rs;
        }

        /*  RQ = 0 (not yet computed) and NQ < lnp. */

        /*  RQ = sqrt(RS) is sufficiently large to (strictly) include NQ nodes. */

        /*  The least squares fit will include NEQ = LNP - 1 equations for */
        /*  5 <= NQ <= NEQ < LMAX <= N-1. */

        if (rq == 0. && nq < lnp)
        {
            neq = lnp - 1;
            rq = sqrt(rs);
            avsq = sum2 / neq;
            /* real ( neq, kind = 8 ) */
        }
        if (nqwmax < lnp)
        {
            goto L4;
        }
        else
        {
            goto L1;
        }

        /*  All LMAX nodes are included in NPTS.   RWS and/or RQ**2 is */
        /*  (arbitrarily) taken to be 10 percent larger than the */
        /*  distance RS to the last node included. */

L3:
        if (rws == 0.)
        {
            rws = rs * 1.1;
        }
        if (rq == 0.)
        {
            neq = lmax;
            rq = sqrt(rs * 1.1);
            avsq = sum2 / neq;
            /* real ( neq, kind = 8 ) */
        }
L4:

        /*  Store RSQ(K), update RSMX if necessary, and compute AV. */

        rsq[k] = rws;
        rsmx = max(rsmx, rws);
        av = sqrt(avsq);

        /*  Set up the augmented regression matrix (transposed) as the */
        /*  columns of B, and zero out the lower triangle (upper */
        /*  triangle of B) with Givens rotations -- QR decomposition */
        /*  with orthogonal matrix Q not stored. */

        i__ = 0;
L5:
        ++i__;
        np = npts[i__ - 1];
        irow = min(i__, 6);
        setup2(&xk, &yk, &fk, &x[np], &y[np], &f[np], &av, &avsq, &rq, &b[
               irow * 6 - 6]);
        if (i__ == 1)
        {
            goto L5;
        }
        i__2 = irow - 1;
        for (j = 1; j <= i__2; ++j)
        {
            jp1 = j + 1;
            givens(&b[j + j * 6 - 7], &b[j + irow * 6 - 7], &c__, &s);
            i__3 = 6 - j;
            rotate(&i__3, &c__, &s, &b[jp1 + j * 6 - 7], &b[jp1 + irow * 6 -
                   7]);
            /*    end do */
            /* L106: */
        }
        if (i__ < neq)
        {
            goto L5;
        }

        /*  Test the system for ill-conditioning. */

        /* Computing MIN */
        d__1 = fabs(b[0]), d__2 = fabs(b[7]), d__1 = min(d__1, d__2), d__2 = fabs(
                                                                                  b[14]), d__1 = min(d__1, d__2), d__2 = fabs(b[21]), d__1 = min(
                                                                                                                                                 d__1, d__2), d__2 = fabs(b[28]);
        dmin__ = min(d__1, d__2);
        if (dtol <= dmin__ * rq)
        {
            goto L13;
        }
        if (neq == lmax)
        {
            goto L10;
        }

        /*  Increase RQ and add another equation to the system to improve conditioning. */
        /*  The number of NPTS elements is also increased if necessary. */

L7:
        rsold = rs;
        ++neq;
        if (neq == lmax)
        {
            goto L9;
        }

        /*   NEQ < LNP. */

        if (neq != lnp)
        {
            np = npts[neq];
            /* Computing 2nd power */
            d__1 = x[np] - xk;
            /* Computing 2nd power */
            d__2 = y[np] - yk;
            rs = d__1 * d__1 + d__2 * d__2;
            if ((rs - rsold) / rs < rtol)
            {
                goto L7;
            }
            rq = sqrt(rs);
            goto L5;
        }

        /*  Add an element to NPTS. */

        ++lnp;
        getnp2(&xk, &yk, &x[1], &y[1], &nnr, &lcell[lcell_offset], &lnext[1],
               &xmn, &ymn, &ddx, &ddy, &np, &rs);
        if (np == 0)
        {
            *ier = 2;
            return;
        }
        npts[lnp - 1] = np;
        if ((rs - rsold) / rs < rtol)
        {
            goto L7;
        }
        rq = sqrt(rs);
        goto L5;
L9:
        rq = sqrt(rs * 1.1);
        goto L5;

        /*  Stabilize the system by damping second partials.  Add multiples of the */
        /*  first three unit vectors to the first three equations. */

L10:
        for (i__ = 1; i__ <= 3; ++i__)
        {
            b[i__ + 29] = sf;
            /*      b(i+1:6,6) = 0.0D+00 */
            for (ijk = i__ + 1; ijk <= 6; ++ijk)
            {
                b[ijk + 29] = 0.;
                /* L400: */
            }
            for (j = i__; j <= 5; ++j)
            {
                jp1 = j + 1;
                givens(&b[j + j * 6 - 7], &b[j + 29], &c__, &s);
                i__2 = 6 - j;
                rotate(&i__2, &c__, &s, &b[jp1 + j * 6 - 7], &b[jp1 + 29]);
                /*      end do */
                /* L108: */
            }
            /*    end do */
            /* L107: */
        }

        /*  Test the stabilized system for ill-conditioning. */

        /* Computing MIN */
        d__1 = fabs(b[0]), d__2 = fabs(b[7]), d__1 = min(d__1, d__2), d__2 = fabs(
                                                                                  b[14]), d__1 = min(d__1, d__2), d__2 = fabs(b[21]), d__1 = min(
                                                                                                                                                 d__1, d__2), d__2 = fabs(b[28]);
        dmin__ = min(d__1, d__2);
        if (dmin__ * rq < dtol)
        {
            *xmin = xmn;
            *ymin = ymn;
            *dx = ddx;
            *dy = ddy;
            *ier = 3;
            return;
        }

        /*  Solve the 5 by 5 triangular system for the coefficients. */

L13:
        for (i__ = 5; i__ >= 1; --i__)
        {
            t = 0.;
            for (j = i__ + 1; j <= 5; ++j)
            {
                t += b[j + i__ * 6 - 7] * a[j + k * 5];
                /*      end do */
                /* L110: */
            }
            a[i__ + k * 5] = (b[i__ * 6 - 1] - t) / b[i__ + i__ * 6 - 7];
            /*    end do */
            /* L109: */
        }

        /*  Scale the coefficients to adjust for the column scaling. */

        for (ijk = 1; ijk <= 3; ++ijk)
        {
            a[ijk + k * 5] /= avsq;
            /* L401: */
        }
        /*    a(1:3,k) = a(1:3,k) / avsq */
        a[k * 5 + 4] /= av;
        a[k * 5 + 5] /= av;

        /*  Unmark K and the elements of NPTS. */

        lnext[k] = -lnext[k];
        i__2 = lnp;
        for (i__ = 1; i__ <= i__2; ++i__)
        {
            np = npts[i__ - 1];
            lnext[np] = -lnext[np];
            /*    end do */
            /* L111: */
        }
        /*  end do */
        /* L105: */
    }

    /*  No errors encountered. */

    *xmin = xmn;
    *ymin = ymn;
    *dx = ddx;
    *dy = ddy;
    *rmax = sqrt(rsmx);
    *ier = 0;
} /* qshep2_ */


/* *****************************************************************************80 */

/* ! QS2VAL evaluates the interpolant function at a point. */

/*  Discussion: */

/*    QS2VAL returns the value Q(PX,PY) where Q is the weighted sum of */
/*    quadratic nodal functions defined by QSHEP2.  If the spatial */
/*    derivatives of Q are also desired, call QS2GRD instead. */

/*    Input parameters are not altered by this function.  The */
/*    parameters other than PX and PY should be input unaltered */
/*    from their values on output from QSHEP2.  This function */
/*    should not be called if a nonzero error flag was returned */
/*    by QSHEP2. */

/*  Modified: */

/*    10 July 1999 */

/*  Author: */

/*    Robert Renka, */
/*    University of North Texas */

/*  Reference: */

/*    Robert Renka, */
/*    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate */
/*    interpolation of scattered data, */
/*    ACM Transactions on Mathematical Software, */
/*    Volume 14, 1988, pages 149-150. */

/*  Parameters: */

/*    Input, double precision PX, PY, the (X,Y) coordinates of the point P at */
/*    which Q is to be evaluated. */

/*    Input, int N, the number of nodes and data values to be */
/*    interpolated.  N must be at least 6. */

/*    Input, double precision X(N), Y(N), the coordinates of the nodes at which */
/*    data has been supplied. */

/*    Input, double precision F(N), the data values at the nodes. */

/*    Input, int NR, the number of rows and columns in the cell grid. */
/*    Refer to subroutine STORE2.  NR must be at least 1. */

/*    Input, int LCELL(NR,NR), the array of nodal indices associated */
/*    with cells.  Refer to STORE2. */

/*    Input, int LNEXT(N), the next-node indices.  Refer to STORE2. */

/*    Input, double precision XMIN, YMIN, DX, DY, the minimum nodal X, Y */
/*    coordinates, and the X, Y dimensions of a cell.  Computed by QSHEP2. */

/*    Input, double precision RMAX, the square root of the largest element */
/*    in RSQ, the maximum radius of influence.  Computed by QSHEP2. */

/*    Input, double precision RSQ(N), the squared radii which enter into the */
/*    weights defining the interpolant Q.  Computed by QSHEP2. */

/*    Input, double precision A(5,N), the coefficients for the nodal functions */
/*    defining the interpolant Q.  Computed by QSHEP2. */

/*    Output, double precision QS2VAL, the interpolated function value */

/*    at (PX,PY). */

double qs2val(double px, double py, int n, double *x,
              double *y, double *f, int nr, int *lcell, int *lnext,
              double xmin, double ymin, double dx, double dy, double rmax, double *rsq, double *a)
{
    /* System generated locals */
    int lcell_dim1, lcell_offset, i__1, i__2;
    double ret_val;

    /* Local variables */
    static int i__, j, k;
    static double w, rd, ds;
    static int kp;
    static double rs, sw, xp, yp;
    static int ijk;
    static double rds, swq, delx, dely;
    static int imax, imin, jmax, jmin;


    /* Parameter adjustments */
    a -= 6;
    --rsq;
    --lnext;
    --f;
    --y;
    --x;
    lcell_dim1 = nr;
    lcell_offset = 1 + lcell_dim1;
    lcell -= lcell_offset;

    /* Function Body */
    ijk = 0;
    xp = px;
    yp = py;
    ret_val = 0.;
    if (n < 6)
    {
        return ret_val;
    }
    else if (nr < 1)
    {
        return ret_val;
    }
    else if (dx <= 0.)
    {
        return ret_val;
    }
    else if (dy <= 0.)
    {
        return ret_val;
    }
    else if (rmax < 0.)
    {
        return ret_val;
    }

    /*  Set imin, imax, jmin, and jmax to cell indices defining */
    /*  the range of the search for nodes whose radii include */
    /*  p.  The cells which must be searched are those intersected */
    /*  by (or contained in) a circle of radius rmax */
    /*  centered at P. */

    imin = (int)((xp - xmin - rmax) / dx) + 1;
    imin = max(imin, 1);
    imax = (int)((xp - xmin + rmax) / dx) + 1;
    imax = min(imax, nr);
    jmin = (int)((yp - ymin - rmax) / dy) + 1;
    jmin = max(jmin, 1);
    jmax = (int)((yp - ymin + rmax) / dy) + 1;
    jmax = min(jmax, nr);

    /*  Test for no cells within the circle of radius RMAX. */

    if (imax < imin || jmax < jmin)
    {
        ret_val = 0.;
        return ret_val;
    }

    /*  Accumulate weight values in SW and weighted nodal function */
    /*  values in swq.  the weights are w(k) = ((r-d)+/(r*d))**2 */
    /*  for r**2 = rsq(k) and d = distance between p and node K. */

    sw = 0.;
    swq = 0.;
    i__1 = jmax;
    for (j = jmin; j <= i__1; ++j)
    {
        i__2 = imax;
        for (i__ = imin; i__ <= i__2; ++i__)
        {
            k = lcell[i__ + j * lcell_dim1];
            if (k != 0)
            {
                /* infinite loop */
                for (ijk = 0; ijk <= 100; ++ijk)
                {
                    delx = xp - x[k];
                    dely = yp - y[k];
                    ds = delx * delx + dely * dely;
                    rs = rsq[k];
                    if (ds < rs)
                    {
                        if (ds == 0.)
                        {
                            ret_val = f[k];
                            return ret_val;
                        }
                        rds = rs * ds;
                        rd = sqrt(rds);
                        w = (rs + ds - rd - rd) / rds;
                        sw += w;
                        swq += w * (f[k] + a[k * 5 + 1] * delx * delx + a[k *
                                    5 + 2] * delx * dely + a[k * 5 + 3] * dely *
                                    dely + a[k * 5 + 4] * delx + a[k * 5 + 5] *
                                    dely);
                    }
                    kp = k;
                    k = lnext[kp];
                    if (k == kp)
                    {
                        /*            exit */
                        goto L214;
                    }
                    /*        end do */
                    /* L114: */
                }
L214:
                ;
            }
            /*    end do */
            /* L113: */
        }
        /*  end do */
        /* L112: */
    }

    /*  SW = 0 if and only if P is not within the radius R(K) for any node K. */

    if (sw == 0.)
    {
        ret_val = 0.;
    }
    else
    {
        ret_val = swq / sw;
    }
    return ret_val;
} /* qs2val_ */


/* *****************************************************************************80 */

/* ! QS2GRD evaluates the interpolant and its first spatial derivatives. */

/*  Discussion: */

/*    QS2GRD computes the value and the gradient at the point (PX,PY) */
/*    of the interpolatory function Q, defined by QSHEP2 for a given set */
/*    of scattered data.  Q(X,Y) is a weighted sum of quadratic */
/*    nodal functions. */

/*    Input parameters are not altered by this subroutine.  The parameters */
/*    other than PX and PY should be input unaltered from their values */
/*    on output from QSHEP2.  This subroutine should not be called if a */
/*    nonzero error flag was returned by QSHEP2. */

/*  Modified: */

/*    10 July 1999 */

/*  Author: */

/*    Robert Renka, */
/*    University of North Texas */

/*  Reference: */

/*    Robert Renka, */
/*    Algorithm 660: QSHEP2D, Quadratic Shepard method for bivariate */
/*    interpolation of scattered data, */
/*    ACM Transactions on Mathematical Software, */
/*    Volume 14, 1988, pages 149-150. */

/*  Parameters: */

/*    Input, double precision PX, PY, the coordinates of the point at which the */
/*    interpolant and its derivatives are to be evaluated. */

/*    Input, int N, the number of nodes and data values which */
/*    are to be interpolated.  N must be at least 6. */

/*    Input, double precision X(N), Y(N), the coordinates of the nodes at which */
/*    data has been supplied. */

/*    Input, double precision F(N), the data values at the nodes. */

/*    Input, int NR, the number of rows and columns in the cell */
/*    grid.  Refer to subroutine STORE2 for details.  NR must be at least 1. */

/*    Input, int LCELL(NR,NR), array of nodal indices associated */
/*    with cells. */

/*    Input, int LNEXT(N), contains next-node indices. */

/*    Input, double precision XMIN, YMIN, DX, DY, the minimum nodal X, Y */
/*    coordinates, and the X, Y dimensions of a cell.  Computed by QSHEP2. */

/*    Input, double precision RMAX, the square root of the largest element */
/*    in RSQ, the maximum radius of influence.  Computed by QSHEP2. */

/*    Input, double precision RSQ(N), the squared radii which enter into the */
/*    weights defining the interpolant Q.  Computed by QSHEP2. */

/*    Input, double precision A(5,N), the coefficients for the nodal functions */
/*    defining the interpolant Q.  Computed by QSHEP2. */

/*    Output, double precision Q, QX, QY, the value of the interpolant, and */
/*    its derivatives with respect to X and Y, at (PX,PY). */

/*    Output, int IER, error indicator. */
/*    0, if no errors were encountered. */
/*    1, if N, NR, DX, DY or RMAX is invalid. */
/*    2, if no errors were encountered but (PX,PY) is not within the */

/*       radius R(K) for any node K and thus Q = QX = QY = 0. */

void qs2grd(double px, double py, int n,
            double *x, double *y, double *f, int nr, int *lcell, int *lnext,
            double xmin, double ymin, double dx, double dy,
            double rmax, double *rsq, double *a,
            double *q, double *qx, double *qy, int *ier)
{
    /* System generated locals */
    int lcell_dim1, lcell_offset, i__1, i__2;

    /* Local variables */
    static int i__, j, k;
    static double t, w, rd, ds;
    static int kp;
    static double qk, rs, sw, xp, yp, wx, wy;
    static int ijk;
    static double rds, qkx, qky, swq, sws, swx, swy, delx, dely;
    static int imax, imin, jmax, jmin;
    static double swqx, swqy;

    /* Parameter adjustments */
    a -= 6;
    --rsq;
    --lnext;
    --f;
    --y;
    --x;
    lcell_dim1 = nr;
    lcell_offset = 1 + lcell_dim1;
    lcell -= lcell_offset;

    /* Function Body */
    xp = px;
    yp = py;
    if (n < 6)
    {
        *ier = 1;
        return;
    }
    else if (nr < 1)
    {
        *ier = 1;
        return;
    }
    else if (dx <= 0.)
    {
        *ier = 1;
        return;
    }
    else if (dy <= 0.)
    {
        *ier = 1;
        return;
    }
    else if (rmax < 0.)
    {
        *ier = 1;
        return;
    }

    /*  Set imin, imax, jmin, and jmax to cell indices defining */
    /*  the range of the search for nodes whose radii include P. */
    /*  The cells which must be searched are those inter- */
    /*  sected by (or contained in) a circle of radius rmax */
    /*  centered at p. */

    imin = (int)((xp - xmin - rmax) / dx) + 1;
    imin = max(imin, 1);
    imax = (int)((xp - xmin + rmax) / dx) + 1;
    imax = min(imax, nr);
    jmin = (int)((yp - ymin - rmax) / dy) + 1;
    jmin = max(jmin, 1);
    jmax = (int)((yp - ymin + rmax) / dy) + 1;
    jmax = min(jmax, nr);

    /*  Test for no cells within the circle of radius RMAX. */

    if (imax < imin || jmax < jmin)
    {
        *q = 0.;
        *qx = 0.;
        *qy = 0.;
        *ier = 2;
        return;
    }

    /*  Q = swq/sw = sum(w(k)*q(k))/sum(w(k)) where the sum is */
    /*  from k = 1 to n, q(k) is the quadratic nodal function, */
    /*  and w(k) = ((r-d)+/(r*d))**2 for radius r(k) and distance d(k).  Thus */

    /*    qx = (swqx*sw - swq*swx)/sw**2  and */
    /*    qy = (swqy*sw - swq*swy)/sw**2 */

    /*  where swqx and swx are partial derivatives with respect */
    /*  to x of swq and sw, respectively.  swqy and swy are */
    /*  defined similarly. */

    sw = 0.;
    swx = 0.;
    swy = 0.;
    swq = 0.;
    swqx = 0.;
    swqy = 0.;

    /*  Outer loop on cells (I,J). */

    i__1 = jmax;
    for (j = jmin; j <= i__1; ++j)
    {
        i__2 = imax;
        for (i__ = imin; i__ <= i__2; ++i__)
        {
            k = lcell[i__ + j * lcell_dim1];

            /*  Inner loop on nodes K. */

            if (k != 0)
            {
                /* infinite loop */
                for (ijk = 0; ijk <= 100; ++ijk)
                {
                    delx = xp - x[k];
                    dely = yp - y[k];
                    ds = delx * delx + dely * dely;
                    rs = rsq[k];
                    if (ds == 0.)
                    {
                        *q = f[k];
                        *qx = a[k * 5 + 4];
                        *qy = a[k * 5 + 5];
                        *ier = 0;
                        return;
                    }
                    if (ds < rs)
                    {
                        rds = rs * ds;
                        rd = sqrt(rds);
                        w = (rs + ds - rd - rd) / rds;
                        t = (rd - rs) * 2. / (ds * rds);
                        wx = delx * t;
                        wy = dely * t;
                        qkx = a[k * 5 + 1] * 2. * delx + a[k * 5 + 2] * dely;
                        qky = a[k * 5 + 2] * delx + a[k * 5 + 3] * 2. * dely;
                        qk = (qkx * delx + qky * dely) / 2.;
                        qkx += a[k * 5 + 4];
                        qky += a[k * 5 + 5];
                        qk = qk + a[k * 5 + 4] * delx + a[k * 5 + 5] * dely +
                            f[k];
                        sw += w;
                        swx += wx;
                        swy += wy;
                        swq += w * qk;
                        swqx = swqx + wx * qk + w * qkx;
                        swqy = swqy + wy * qk + w * qky;
                    }
                    kp = k;
                    k = lnext[kp];
                    if (k == kp)
                    {
                        /* exit */
                        goto L204;
                    }
                    /*        end do */
                    /* L104: */
                }
L204:
                ;
            }
            /*    end do */
            /* L103: */
        }
        /*  end do */
        /* L102: */
    }

    /*  SW = 0 if and only if P is not within the radius R(K) for any node K. */

    if (sw != 0.)
    {
        *q = swq / sw;
        sws = sw * sw;
        *qx = (swqx * sw - swq * swx) / sws;
        *qy = (swqy * sw - swq * swy) / sws;
        *ier = 0;
    }
    else
    {
        *q = 0.;
        *qx = 0.;
        *qy = 0.;
        *ier = 2;
    }
} /* qs2grd_ */


